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辛 算法 作为 研究 哈密 顿 系统 长 期 定性 演化 的 最 佳 积 ， 自 问世 以 来 就 受到 了 很 大 的 


关注 。 通 过 对 哈密 顿 函数 的 截断 误差 分 析 ， 可 以 从 不 同 角度 构造 出 较 高 精度 的 辛 算 法 ， 也 可 以 通 


过 引入 


正规 化 技术 实现 自动 调整 积分 步 长 和 改善 数值 稳定 性 。 从 辛 算法 的 表现 形式 可 以 将 它 分 


AY SENURI ES P 


VER. 当 哈 密 顿 系统 能 够 分 解 为 儿 个 可 积 部 分 且 每 部 分 的 解 能 用 时 间 显 函数 来 表 


示 时 ， 可 以 构造 显 式 算 法 。 显 式 算法 有 非 力 梯度 显 式 辛 算法 、 力 梯度 辛 算法 、 辛 校正 、 类 高 阶 辛 
算法 四 种 。 当 哈密 顿 系统 变量 不 能 分 离 时 ， 适 合 应 用 隐 式 辛 算法 和 扩充 相 空间 对 称 算 法 求解 。 分 
别 对 这 些 算 法 的 构造 方法 及 其 适用 的 物理 模型 进行 归纳 对 比 ， 分 析 了 各 种 辛 算 法 的 优 劣 性 和 发 


Xx 键 
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哈密 顿 系统 从 物理 本 质 上 具有 辛 结构 的 不 变性 。20 世纪 80 年 代 初 期 ， 中 国学 者 冯 康 先 
ET 与 Ruth” 几 乎 同时 提出 能 够 保持 哈密 顿 相 流 辛 结构 的 数值 积分 算法 ， 解 决 了 传统 算法 
因 长 期 对 时 间 积分 不 能 保持 系统 的 能 量 守恒 的 问题 "。 前 者 主要 对 不 可 分 的 哈密 顿 系统 采用 


H 


展 趋势 ， 


对 如 何 选择 辛 算法 高 效 高 精度 地 解决 实际 问题 提供 了 一 定 的 理论 和 数值 计算 依据 。 


词 : PEA; PAWR; FRE; JARRA: PRA AR RIE 
中 图 分 类 号 : P138 文献 标识 码 : A 


A 


以 隐 式 中 点 法 为 基础 来 构建 高 阶 隐 式 辛 算法 ， 后 者 主要 是 针对 可 分 解 为 动能 了 和 势能 站 的 
哈密 顿 系统 建立 显 式 辛 算法。 显 式 辛 算法 和 隐 式 辛 算法 的 区 别 在 于 积分 过 程 中 是 否 需 要 迭 


代 计算 。 


在 Ruth 建立 的 开 十 六 形式 的 显 式 辛 算法 的 基础 上 ， 高 阶 辛 算法 的 研究 和 应 用 得 到 了 
快速 发 展 ”。 当 哈密 顿 分 解 成 主要 的 未 受 摄 部 分 Ho 和 次 要 的 摄 动 部 分 eH, 且 两 部 分 均 可 积 
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时 ”， 可 以 通过 升 高 摄 动 项 eH, 的 小 参数 = 的 阶 提高 精度 ， 即 类 高 阶 辛 算法 (pseudo-high- 
order-symplectic integrator, PSD” 和 辛 校正 (wisdom-holman-touma correction, WHT)". 
力 梯度 辛 算法 "是 在 算法 的 Lie 算 子 中 嵌入 了 力 梯 度 算 子 ， 可 以 有 效 避 免 负 时 间 系 数 的 
HI Ruth” f Chin™ ”分 别 构造 了 三 阶 力 梯度 辛 算法 和 四 阶 力 梯度 辛 算法 。 李 荣 和 伍 
ATT gA McLachlan ™ 的 算法 优化 思想 构造 了 三 阶 优化 力 梯度 辛 算法 和 四 阶 优化 力 梯度 
辛 算法 。 陈 云龙 和 伍 坎 ”论证 了 力 梯度 辛 算法 在 求解 旋转 坐标 系 下 的 圆 型 限制 性 三 体 问题 
的 有 效 性 。 
哈密 顿 系 统 变量 不 可 分 离 时 ， 显 式 辛 算法 一 般 不 能 直接 应 用 ， 隐 式 辛 算法 应 该 适合 应 
j. MAMAS RIES Kee eth ET T. Be eee i 
系统 ， 但 计算 效率 低 ， 不 宜 作 为 首选 数值 方法 。 为 使 显 式 辛 算法 在 不 可 分 哈密 顿 系 统 中 得 到 
应 用 ，Pihajoki 在 不 可 分 哈密 顿 系统 中 增加 了 一 组 与 原 相 空间 动量 坐标 相同 的 量 ， 提 出 了 
相 空 间 扩 充 对 称 算法 ， 但 系统 的 辛 结构 被 破坏 。Tao” 在 Pihajoki 的 基础 上 将 哈密 顿 系 统 分 
成 三 部 分 ， 第 三 部 分 是 两 组 坐标 动量 差 的 平方 组 成 的 约束 哈密 顿 ， 保 证 了 数值 解 在 扩展 相 空 
间 中 的 辛 结构 。 在 此 基础 上 ，Liu 和 Wu” 构造 了 连续 坐标 动量 置换 相 空间 扩大 法 ， 算 法 
适合 求解 后 牛顿 哈密 顿 问题 和 旋转 哈密 顿 问 题 等 哈密 顿 系统 不 可 分 的 问题 ， 但 在 旋转 致密 
双星 后 牛顿 哈密 顿 系统 中 不 能 保持 系统 的 能 量 稳定 。Luo 等 人 ”提出 了 中 点 置换 相 空间 
扩大 方法 ， 这 种 方法 优越 于 Pihajoki 方法 和 刘磊 等 人 的 方法 ， 它 对 算法 的 算 子 个 数 没 有 要 
求 ， 适 用 范围 更 广 。Wmu 等 人 ”提出 的 扩大 相 空间 优化 Forest-Ruth 算法 在 平面 圆 型 限制 性 
三 体 问题 和 致密 双星 问题 中 与 未 优化 算法 相 比 精度 得 到 了 明显 改善 ， 且 表明 采用 中 点 置换 
能 使 优化 的 Forest-Ruth 算法 具有 更 好 的 数值 表现 。Li 和 Wu” 在 Mikkola 等 人 ”的 基础 
上 结合 扩大 相 空间 中 的 几 种 置换 方法 ， 提 出 了 扩大 相 空 间 的 对 数 哈密 顿 显 式 对 称 法 。 这 些 方 
去 在 圆 型 限制 性 三 体 问 题 、 三 阶 后 牛顿 自 旋 致密 双星 系统 和 Ernst-Schwarzschild 黑洞 等 模 
型 中 表现 出 了 高 精度 和 高 效率 。 辛 算法 在 解决 各 类 天 文学 物理 模型 中 得 到 了 较为 广泛 的 应 
o 本文 的 主要 目的 是 对 上 述 辛 算法 及 类 辛 算法 的 发 展 、 构 建 以 及 适用 物理 模型 进行 阐述 。 


Re 


^H 


2 六 算法 的 建立 
按照 工 十 V 分 解 形式 ， 假 设 哈密 顿 可 以 表示 为 ; 


Hp, d - T) V() = BE +v), () 


式 中 ，p, q 分 别 表示 广义 坐标 和 广义 动量 。 相 应 的 正则 运动 方程 : 


d= = (2 1,2,8,2-en) (2) 


或 取 z = (p, q) 
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Koh, Ly Æ Li RT, PHD AM BATT AV 的 Lie 导数 : 


n o 
A=4{,T}=} pi 
a” oqi , (4) 


n o 


Ks fi = 一 站 是 函数 的 第 i 个 分 量 。 分 别 将 运动 方程 的 解 441 和 popi 在 dn 处 泰勒 展 
开 ， 则 辛 算法 可 以 表示 为 动能 和 势能 所 对 应 的 Lie 算 子 的 指数 积 所 构成 的 一 系列 组 合 : 


| dui = e (AtB)g, 


paries er (^*B)p,, , (5) 


AF, 7 为 时 间 步 长 。 设 W = {, H} 将 W WAR MHA AGI BRS RAG T 
成 ， 即 
ors II etter O+ (pH) um Or t+) . (6) 


2 


a; 和 访 是 由 特定 的 阶 确定 的 系数 ，O(A<+1) 指 哈密 顿 函数 的 IK 阶 截断 误差 。 若 使 Lie HF 
的 个 数 等 于 阶 条 件 的 个 数 ， 时 间 系 数 可 以 很 容易 被 算出 。Ruth 在 文献 中 提出 了 比较 常用 
的 二 阶 和 三 阶 辛 算法 。 在 此 基础 上 ，Forest-Ruth 四 阶 非 力 梯 度 辛 算法 ”和 Yoshida 高 阶 非 
力 梯度 辛 算法 ”也 相继 被 建立 。 


3 AEA BRE Se SUE SUE 


3.1 T+V 分 解法 
3.1.1 未 优化 的 非 力 梯 度 显 式 辛 算法 
按照 Ruth 的 方法 ， 对 于 可 积分 离 的 哈密 顿 系统 ， 辛 算法 构造 如 下 ， 
(1) 一 阶 显 式 辛 算法 


(2) 二 阶 显 式 对 称 辛 算法 
二 阶 显 式 对 称 辛 算法 由 三 个 单 指数 算 子 构成 ， 它 的 构造 形式 为 : 


| Sox : e37AqQT Be 37A = e'W 


Soy : est BeTAgazTB = eTW 


(3) 三 阶 显 式 辛 算法 
三 阶 显 式 辛 算法 由 6 个 单 指数 算 子 构成 ， 它 的 构造 形式 为 S3: 


A.B. _27r4 37 274 27B 
er4e 到 7Be 574ei7Bes7A4e2 =—e7W . (9) 
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(4) 四 阶 显 式 Forest-Ruth 辛 算法 (fourth-order-Forest-Ruth symplectic integrator, FR) 
Forest-Ruth 算法 从 式 (6) 出 发 ， 它 由 了 个 单 指数 构成 S4: 


e47c/2eBrce4(I-c)7/2eB(I-2c)re4U-c)r/2eBrce4rc/2 ，C 二 1/(2 v2) . (10) 


构造 高 阶 算法 的 另 一 个 思路 是 直接 对 二 阶 辛 算法 进行 多 重 积 ，Yoshida” 在 Ruth" 的 理 
论 基 础 上 用 三 个 二 阶 辛 算法 对 称 组 合 得 到 四 阶 显 式 Yoshida 辛 算法 Syo: 


Syo = So(cT)S2((1 — 2c)7)S2(cT) . (11) 


同年 Forest 和 Ruth 提出 了 Forest-Ruth 四 阶 辛 算法 ， 这 两 种 辛 算法 构造 可 以 达到 等 
价 的 效果 。 按 照 上 述 算法 构造 方式 ， 高 阶 辛 算法 95,19 可 以 被 构造 ”， 


Son+2 = Son(aT)Son((1 — 2a)7)S2n(aT), a=1/(2 -2/21D)) , (12) 


3.1.2 ”优化 的 非 力 梯度 显 式 辛 算法 

McLachlan 1992 年 提出 了 算法 的 优化 概念 ， 通 过 将 主要 截断 误差 项 的 影响 降 到 最 小 ， 
得 到 优化 系数 ， 从 而 构造 优化 的 三 阶 显 式 辛 算法 和 三 阶 显 式 辛 算法 。 在 此 基础 上 ， 四 阶 优 化 
非 力 梯度 算法 等 高 阶 优化 算法 也 被 建立 起 来 。 

(1) 二 阶 优化 的 显 式 辛 算法 和 三 阶 优化 显 式 辛 算法 
单一 化 系数 表示 的 二 阶 显 式 辛 算法 的 主要 误差 函数 为 : 


1 — 4a3 


1 
ho = F'P?-( | Fr? p’ 1 
2 2 3a3) 8(1 — a2) , ( 3) 
E aV OT oua - "t 
规定 让 = 一己 , P= > F' 和 P 表示 两 者 的 导数 ，as 表示 截断 误差 项 系数 。McLachlan 
q Dp 


4 ha 的 系数 平方 和 最 小 ， 得 到 两 个 极 小 值 ag = 1 e a» 的 值 越 小 误差 就 越 小 。 与 未 
优化 的 二 阶 显 式 辛 算法 和 类 高 阶 辛 算法 的 截断 误差 常数 相 比 ， 二 阶 最 优化 算法 精度 提高 了 
0.39， 如 表 1 所 示 。 

HE, McLachlan 将 ?” 阶 辛 算法 的 主要 误差 函数 表示 为 ; 


-ÈA 2 (a, b)gi(F, P)|¢ao.p0) d (14) 


f; 和 gi 表示 误差 常数 。 令 该 式 的 系数 平方 和 最 小 ， 得 到 三 阶 最 优化 辛 算法 Os: 


g^17AgbiTBga2TAgbar Bast AgbsT B . (15) 


(2) 四 阶 显 式 优 化 Forest-Ruth 算法 (fourth-order-optimal-Forest-Ruth integrator, OFR) 
Omelyan 4& A" 将 四 阶 显 式 Forest-Ruth 算法 拓展 分 解 为 : 


SA BEPe4(1-2AX)Pn/2eBXxhe4XheB(I-2(X+T6))ne4AXheBXhe4(I-2AX)Pn/2eBSP | (16) 
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Xi 几 种 显 式 辛 算法 的 截断 误差 比较 
oe a b 误差 常数 
二 阶 显 式 辛 算法 al = a» = 1/2 b; =0, b2 =1 0.280 (0.070) 
类 二 阶 辛 算 法 ai =1, a» =0 b; = b2 = 1/2 0.280 
、 E EM 2 2 
优化 的 二 阶 显 式 辛 算法 a= lila bo Y ilb 0.043 
_. 、 2 2 7 3 1 
— WMV RH EA ya 5 
三 阶 显 式 辛 算法 a1 = 3,02 3 03 1 bi 278 b2 a bs m 0.113 
= 2 1 
优化 的 三 阶 显 式 辛 算法 dm. cA 
ia (7 RE 
w= 26 1 +z 
3 9z 
y= 70 +u’) 
(i _wv REN 
Bre g + Vy) 3/U 
—0.919 661 523 017 399 857 
1 ay, 1 ay 本 | 
a2 dan 2 5 a2 Aai 2 bi = QA4—i 0.058 


Cy EX, 入 为 3 个 步 长 系数 ， 它 们 之 间 存 在 两 个 约束 条 件 ，4 的 系数 alx A) = 0，B 的 系数 
= BEX A) = 0。 由 于 仍 存在 一 个 自由 度 ， 步 长 系数 的 取 值 将 有 多 种 解 。 从 McLachlan 的 优 
Te 化 思想 出 发 ， 令 四 阶 显 式 Forest-Ruth 算法 的 三 阶 截 断 误 差 为 零 及 五 阶 截断 误差 项 系数 平方 
和 最 小 ， 得 到 最 优 解 : 


YD £ = 0.172 086 559 029 514 3 , 
c A = —0.091 562 030 755 156 78 , 
C x = —0.161 621 762 210 7222 . 
= 该 组 合 系数 解 下 的 OFR 算法 的 截断 误差 值 为 yain = 0.000 92， 四 阶 FR 辛 算法 的 截断 误差 
值 为 Yea = 0.039， 优 化 后 的 精度 得 到 了 显著 提高 
c (3) 优化 的 四 阶 类 Suziki 算法 (fourth-order-pseudo-Suziki integrator, SU) 

Omelyan 等 人 ”在 Yoshida 构建 的 四 阶 辛 算法 两 端 各 加 入 两 个 相同 算 子 构建 了 SU A 
法 ， 它 由 5 个 算 子 组 合 而 成 : 


Ssu = S»(aih)S»(ash)S»((1 = 2a, = 2a2)h)S2(a2h)S2(arh) ; (17) 
a, = 0.322 137 596 081 798 4, aa = 0.541 316 548 170 043 0. 


考虑 满足 阶 条 件 和 截断 误差 h^ 系数 平方 和 最 小 ， 得 到 的 哈密 顿 截断 误差 值 为 y = 0.001 1, 
与 OFR 算法 相 比 ， 误 差 精度 低 一 个 量 级 。 
3.2 H = Hy eH, 分 解法 

H+, Ho MA, AAR, Ho 是 主要 项 ，Hi 是 摄 动 项 ， 小 参数 s 表示 HQ 与 Ho AEK e 
的 差别 ， 具 体 分 解 过 程 在 文献 [35] 中 有 详细 描述 。 
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规定 4 = {, Ho} MB = {,My}, 由 辛 算法 的 构造 公式 erW = erAerB" 并 反复 利用 
Baker-Campbell-Hausdroff (BCH) 公式 ”得 到 哈密 顿 函数 的 误差 表达 式 : 


Hae = e(Q1T 二 + )+E (DT eer) s (18) 


ay, by 是 误差 项 系数 ， 从 上 式 可 以 得 到 三 种 提高 辛 积 分 器 的 精度 方法 : (1) 提高 7 A RK 
得 到 高 阶 辛 算法 ，(2) 改进 哈密 顿 的 分 解 方式 和 构造 方式 ，s 越 小 ， 精 度 越 高 ， 比 如 哈密 顿 
函数 握 动 分 解 "、 时 间 变 换 以 及 扩展 相 空 间 ; (3) 升 高 = 的 寒 次 ， 得 到 辛 校正 ”。 若 误差 项 
中 的 O(e2r"+1) 仍 存在 ， 升 高 r 的 宕 次 得 到 类 高 阶 辛 方法 ”。 这 两 种 类 辛 算法 不 需要 减 小 
积分 步 长 ， 也 不 用 升 高 积分 的 阶 就 能 优化 积分 精度 和 效率 。 一 般 情 况 下 ， 哈 密 顿 函数 按照 
H = Ho + eH, aN, ST+V 分 解 相 比 ， 算 法 精度 有 显著 提高 。 如 图 1 所 示 ， 以 纯 开 普 
勒 问题 为 例 ， 摄 动 分 解 时 的 算法 的 数值 精度 比 了 十 V 分 解 提 高 了 接近 5 个 量 级 。 


lg(AE/E) 
( 


lgt 
0.999 0.001 mx 


T 


注 ， 二 体 问题 摄 动 分 解 形式 是 H p? 


数值 精度 提高 了 接近 5 个 量 级 。 


1 ， 纯 开 普 勒 问题 在 两 种 哈密 顿 分 解 形式 中 的 能 量 误差 


， 步 长 取 周 期 的 1/90, 5 T +V 分 解 形式 相 比 ， 


3.21 d EE 
Wisdom 和 Holman™ 在 雅 可 比 坐 标 系 建立 了 如 下 的 二 阶 辛 算法 (second-order-wisdom- 
verlet, MV2): 


A+B = ehNB -— et B obi B gao AGbo B . (19) 


AP, N =1, bo = bi = 1/2, ay = 1. MV2 的 数值 结果 对 应 的 哈密 顿 函 数 太 和 真实 哈密 顿 
函数 H 的 关系 为 所 = H+ Has Wisdom 等 人 "运用 Lie 级 数 算 子 构造 一 个 正则 变换 的 生 
成 函数 Ws Hor 不 含 摄 动 小 参数 。 的 一 或 者 二 次 寡 项 ， 从 而 构造 辛 校正 公式 。 辛 校正 具有 
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如 下 特点 : (1) 在 不 需要 减少 步 长 也 不 需要 升 高 阶 数 的 情况 下 ， 通 过 消去 误差 项 中 的 = 或 e? 
项 ， 实 现 一 阶 或 者 二 阶 辛 校正 来 显著 提高 数值 精度 ; (2) 辛 校正 不 严格 保 辛 ， 但 没有 引入 耗 
散 机 制 ，(3) 计算 效率 高 。 

WHT 计算 生成 函数 W 的 过 程 较为 繁琐 ， 而 且 校 正 子 C 只 适用 于 上 述 辛 算法 ， 不 具备 
一 般 性 。Mikkola 和 Palmer" " 借助 Euler-Mclaurin 式 推导 出 辛 校正 的 生成 函数 W. A, 
Chambers 和 Murisun" 建立 了 类 高 阶 辛 算法 。Duncan 等 人 号 号 借助 日 心 坐标 系 将 太阳 系 
动力 学 中 的 哈密 顿 函 数 分 为 可 积 的 三 部 分 。Malhotra 加 研究 伽利略 卫星 的 潮汐 演化 时 ， 将 
哈密 顿 系统 分 为 了 14 个 可 积 部 分 。Wu 等 人 号 号 令 哈 密 顿 函数 的 多 部 分 分 解 更 加 一 般 化 ， 
将 辛 算 法 中 的 名 分 成 了 NN 个 子 步 ， 并 用 数值 方法 讨论 了 哈密 顿 在 具有 N41 个 可 积 部 分 时 
的 算法 精度 (主要 是 一 阶 和 二 阶 )。 

Wu 等 人 ™ 利用 伯 努 利多 项 式 ”推出 辛 算法 的 哈密 顿 的 误差 函数 : 从 哈密 顿 函数 误差 项 
HR, Wee Me? 项 ， 直 接 确定 校正 子 C 而 不 计算 生成 函数 WV， 分 别 得 到 一 阶 辛 校正 和 二 
阶 辛 校正 。 

(1) e 的 一 阶 辛 校正 

由 误差 函数 得 到 一 阶 辛 校正 算 子 ” ， 从 而 消去 截断 误差 项 中 e; 的 一 阶 项 ， 一 阶 校正 子 
的 表达 式 为 : 


1 Be) (1) 2j—1 


Cs, =7 (Bi (2 J T 9271 B;) 
i=1 j=1 J): 
C $ 3 a BCD (0) 2jg2j-1p (20) 
== €iT i * 
Sox A (25)!192-1 vA 
N n BO). Sous 

C == - e,T I S) B; 

Soy »»» (25)! A 


Hit, A= Lm» Bi = Lyp H% F [A,B] = AB — BA, [A,B,C] = [A, [B,C]; Mid 
[A, B] =S4B, [A, A, B, C,C] 2 939 5S6, B(x) 是 伯 努 利 公 式 ， 满 足 递 推 公式 : 


m O (x gm 
2 a pens en 


i= 


式 (20) 适用 于 有 n 个 可 积 部 分 的 哈密 顿 系统 的 任意 阶 辛 算 法 ， 哈 密 顿 的 截断 误差 可 以 通过 
matlab 或 mathematica 等 数学 软件 直接 确定 ”， 提 高 了 计算 效率 。 

(2) e 的 二 阶 辛 校正 

Wu 等 人 讨论 了 二 阶 显 式 辛 算法 Sox 的 二 阶 辛 校正 过 程 ， 该 推导 过 程 适用 于 一 般 n 
阶 辛 算法 。.52x 关于 e 的 校正 公式 如 下 : 


e —C ,— C» K CoC 
Sox =e "e "?el ee" , (22) 
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其 中 ， 
/ T I T 
ek = FER ez4 ; 
N N N i A i 
1 2 3 
Be 2 &Bi + 2 p» (NEiEIT Sp,SaAB; — 14406i€3T Sp, S4 D; Tee :) ; 
t= i=1 j=1 
3 i) 1 i i +1 
GS cilen TY Bi t chya T S UNIS ug Bite) te, 
i=1 
_ N N i i 
z = x C2 
C» £z pcr , (5760137 S 5,94 B; +) » 
I= J= 


C; 可 用 @ 近似 代替 ” 。 在 太阳 -木星 -土星 模型 中 ， 从 能 量 误差 、 计 算 效 率 和 平均 经 度 误 差 
方面 进行 数值 比较 ， 辛 校正 的 计算 精度 比 显 式 辛 算法 提高 了 一 个 量 级 并 具有 更 好 的 数值 稳 
定性 。 
3.2.2 ”类 高 阶 辛 算法 

类 高 阶 辛 算法 是 在 摄 动 分 解 的 基础 上 形成 高 效 的 积分 方法 。 相 比较 高 阶 辛 算法 的 子 步 
随时 间 步 长 的 增加 而 迅速 增加 ， 类 高 阶 辛 算法 的 子 步 更 少 且 计 算 效 率 更 高 。Chambers 和 
Murison" 从 四 阶 和 六 阶 辛 算法 出 发 ， 通 过 构造 合适 系数 使 误差 项 中 不 含 ert 或 sr6， 构 造 过 
FEU F: 


e22 7Aebl eTBeai 7Aebl <sTBea2 TA 


@(a1+2a2)TA+2bie7B H ec? (n )[(a14-2a2)? —6a2 (a1--a2)]S92 B 


23 
e?r CE) a a1)82A er5 (PL )[(a4 -2a5)* 3022 (ai Fa3)?]S2 B ( ) 
e475 (PE) (1600 7a1)S4A+-- 
当 
a, +2a2=1, 26-1, 1-—6ag(l—az)=0. 
时 ， 可 消去 eT? 项 。 
` Ny B DEN Me s+ | 
文中 列举 了 类 四 阶 辛 算法 和 类 六 阶 辛 算法 “: 
PSIA = e$ 40- We E Bea AE Bet A) 
— grae (Ss) A- ds 95 BH l (24) 
PSI6 = e$40- Ya) eS Bose Ao Bose AQ te Betal- Je) 
— erue URS A OCT) (25) 


八 阶 和 十 阶 等 高 阶 算法 也 可 被 构造 。 从 算法 的 构成 特点 看 ， 类 高 阶 辛 算法 能 保持 系统 能 量 稳 
定 下 ， 计 算 效率 也 明显 优 于 同 阶 的 通常 辛 算法 ™“”。 从 算法 的 截断 误差 项 来 看 ， 它 可 以 
看 作 是 对 二 阶 辛 算法 的 部 分 校正 。 类 高 阶 辛 算法 到 达 一 定 阶 时 ， 精 度 将 不 再 随 阶 数 的 增加 而 
提高 。Laskar 和 Robutel” 指出 类 六 阶 或 八 阶 的 数值 效果 最 好 。Wnu 等 人 ”在 日 - 木 - 土 三 体 
模型 中 对 类 四 阶 辛 算法 PSI4、 它 的 校正 算法 及 通常 四 阶 辛 算法 54 三 种 算法 进行 数值 模拟 ， 
结果 表明 三 者 的 能 量 误差 精度 表现 为 同一 量 级 。 
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4 Fai BEE SL 


通常 辛 算法 在 积分 过 程 中 难免 会 出 现 负 积分 步 长 ， 但 在 不 可 逆 的 力学 问题 中 需 使 得 积 
分 步 长 全 为 正 。Ruth” 在 三 阶 辛 算法 的 算 子 组 合 中 引入 力 梯度 算 子 [V, P V], aM 
分 子 步 的 系数 均 为 正 。Chin 等 人 ”把 力 梯度 辛 算法 发 展 到 四 阶 ， 与 同 阶 的 非 力 梯度 
辛 算法 相 比 精度 提高 了 10 ~ 80 倍 ”， 在 圆 型 限制 性 三 体 问 题 、 含 时 引力 场 问题 ”及 量子 
力学 问题 “中 得 到 了 很 好 的 应 用 。Sun & AP" E Chin 提出 的 四 阶 力 梯度 辛 算法 基础 上 构造 
了 两 种 含有 三 阶 导数 项 的 四 阶 力 梯 度 辛 算法 ， 在 Hénon-Heiles 系统 、 四 极 矩 核 - 壳 模 型 以 及 
含 扁 率 限制 性 三 体 问题 中 进行 数值 模拟 ， 精 度 和 计算 效率 都 得 到 了 提高 。Xu 和 Wi wy 
] 这 种 力 梯度 辛 算法 求解 摄 动 二 体 问 题 和 雅 可 比 坐 标 系 下 的 摄 动 N 体 问题 。Omelyan 等 
人 “四 将 最 优化 思想 引入 到 力 梯度 辛 算法 ， 按 照 对 称 组 合 的 形式 构造 了 三 阶 和 四 阶 最 优 
化 力 梯度 辛 算 法 ， 并 将 这 两 种 方法 运用 到 天 体力 学 、 分 子 动力 学 和 量子 力学 等 领域 ， 与 通 
常 辛 算法 相 比 ， 数 值 精度 得 到 了 明显 的 改善 ， 但 是 步 长 系数 还 是 不 可 避免 地 用 到 了 负数 。 
李 荣 等 人 也 中 构造 了 步 长 系数 全 为 正 步 长 的 不 对 称 的 三 阶 、 四 阶 最 优化 力 梯度 辛 算 
法 ， 这 种 方法 在 求解 摄 动 开 普 勒 混沌 问题 ” 的 能 量 精度 和 定 态 薛 定 雇 方程 能 量 本 征 值 问题 
ri 7207 ae gp ag gi pr RECS. Wes JERV IC 论证 了 质心 旋转 坐标 系 中 动能 全 不 是 动 
E 己 的 严格 二 次 型 时 ， 力 梯度 辛 算法 仍然 适合 求解 。 


Seen 


定义 力 梯度 算 子 为 : 
3 a 
[C] = [B, A, B] = {, {{V,T}, V} = PON 
asa ey vie (26) 
ij-l '0q; Opi i=1 Uo Opi" 
OV OV D29 
其 设 = a a 将 EA. 
组 合 可 构造 出 力 梯 度 辛 算法 : 
e W - eT (4t B) = eair4ebirB+cirC . (27) 
A4. 非 优化 力 梯度 辛 算法 
(1) 三 阶 力 梯 度 辛 算法 F3 
egr4esirBesr4eirB+ 柱 rC =o | (28) 


这 种 算法 的 积分 步 长 都 为 正 ， 但 力 梯 度 算 子 C 的 计算 过 程 较 复杂 ， 早 期 很 少 受到 重视 。 
(2) 四 阶 力 梯度 辛 算法 F4 


Chinay iN AdE 
ChinaA IVE 1 F HH Fl 
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Chin 及 其 合作 者 在 Ruth KÆ E35 Z1 b6 RE SE SEA Js gp pup n mem, 


A: e$ 7Be3TAgo TB iT Cog TAgaTB Z e'W (29) 
lg 3 Es 3 afjl 
B: e20 FTA GST B+ dg (2— V3)r Ce 374937 Bt ag (2- V3)r Cea FTA = ew (30) 
1 3 i i Ql 73 i 3 1 
C: es Aes BegtTAgaTBt 937 Ces74eg7Be674 = e'W (31) 
3 3 
D: eat Bt saat Ce3TAgs Bg T Ags TB Q3 AgaT Bt aaiT C... e'W (32) 


与 通常 四 阶 辛 算法 相 比 ， 精 度 提 高 了 10 80 fi. 
Xu 和 Wu” 按照 Yoshida 构造 标准 高 阶 辛 算法 的 思路 将 算 子 C 与 A, B 按照 两 种 不 同 
方式 对 称 组 合 得 到 了 两 类 新 的 四 阶 力 梯度 辛 算法 。 
(a) 力 梯度 算 子 蔡 套 在 正中 间 时 ， 表 示 为 : 


< 


3 
eal74epl7Bea274eba7 刀 bs7 Ceaa7r4epl7Beal74 = ew (33) 


ew (b) FIBRES KBE Bam, KRN: 


3 . 3 
ghirB bar C ea1T AgbaT B ga27 AgbaT B gayT Agbyr Bt bet C = eV . (34) 


根据 BCH 公式 ， 从 上 式 组 合算 子 的 中 间 开 始 实施 循环 推导 得 到 W 的 具体 表达 式 ， 由 
最 优化 思想 得 到 一 系列 新 的 正 积分 系数 ， 如 表 2 和 表 3 所 示 。 A 和 By 型 算法 对 应 于 Chin 
的 四 阶 C 和 D 算法 。Xu 和 Wu” 用 这 两 种 新 的 四 阶 力 梯度 辛 算法 求解 摄 动 开 普 勒 问题 和 雅 
可 比 坐标 系 下 的 N 体 问题 ， 以 及 深入 研究 了 行星 磁 气 圈 内 的 带电 粒子 的 混沌 运动 。 如 图 加 
所 示 ， 两 种 分 解 形式 中 力 梯度 辛 算法 的 数值 性 能 没有 明显 差异 ， 摄 动 分 解 形 式 中 算法 精度 有 
显著 提高 。 摄 动 分 解 形 式 中 的 力 梯度 算 子 [Ho Hi, Ho] 的 计算 比 计算 Ho 更 加 简单 ， 与 Ha 
相 比 不 会 增加 大 量 额外 的 计算 成 本 ， 所 以 是 非常 值得 推荐 其 用 于 研究 雅 可 比 坐 标 系 下 的 N 


体 哈密 顿 问题 一 ” 。 


#2 A 型 四 阶 力 梯度 辛 算法 的 几 组 正 积分 系数 


= 算法 A» As A4 
1 V15 1 V2 
mS eoe .181 441 601 ft eS 
a1 3 TE 0.18 601 770 87 hw. 
az vis 0.318 558 398 229 129 v2 
12 4 
b 2 0.410 592 148 47 E 
J 5 i 3 
2 1 
be - 0.178 815 703 059 189 B 
1 V15 1 v2 
bs 2c 0.006 240 214 404 679 3 Boe 
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表 3 ”B 型 四 阶 力 梯度 辛 算法 的 几 组 正 积分 系数 
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算法 B» Ba Ba 
ay 3 0.399 986 828 125 39 0.409 715 409 973 941 
Q2 : 0.200 026 350 374 923 0.1 805 691 800 511 07 
bi 元 0.152 773 965 219 889 0.155 431 946 448 732 
b2 zu 0.003 279 056 273 196 9 0.003 488 836 809 494 1 
b3 = 0.347 226 034 780 111 0.344 568 053 551 268 
X10? /  -( /—» X105 
e-FR 
os 1.0 
0.4 E 
RI 
0.0 <q 0.0 
-0.4 
-1.0 
0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 


t/104 t/104 
a) c) 
X105 X10? 
0.2 eAd -0.2 5-B2 
-0.6 
= -0.2 B 
< a -1.0 
06 -1.4 
0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 0.0 0.2 04 0.6 0.8 
t/104 t/10! t/10* 
d) e) f) 
VE: 图 中 ce 一 表示 摄 动 分 解 ， 未 标注 表示 下 十 V Dii ET +V 型 哈密 顿 分 解 中 ， 新 的 四 阶 力 梯度 辛 算法 的 精 
度 远 远 高 于 标准 四 阶 辛 算法 ， 哈 密 顿 摄 动 分 解 中 ， 就 平 经 度 和 相对 位 置 精 度 而 言 几乎 与 标准 四 阶 辛 算法 相同 。 


图 2 FRIA, B 两 种 类 型 四 阶 力 梯度 辛 算法 在 两 种 哈密 顿 分 解 方 式 中 的 能 量 误差 


4.2 ”优化 力 梯度 辛 算法 
(1) 三 阶 最 优化 力 梯度 辛 算法 


Lie 算 子 系数 的 求解 只 适用 于 对 称 辛 算法 ， 给 不 对 称 的 力 梯 度 辛 算法 的 构造 增加 了 
难 。 李 荣 等 人 构造 了 不 对 称 的 三 阶 和 四 阶 最 优 力 梯度 辛 算法 。 根 据 Chin 的 思路 ， 
中 构造 了 如 下 的 不 对 称 算 子 组 合 : 


e 


根据 BCH 公式 得 到 阶 条 件 的 关系 式 ， 由 最 优化 思想 ， 即 加 入 附加 的 约束 条 件 ” ( 令 其 四 


s 
文 


TC 
bit B+k eal74epl7Bea274eb27 刀 = e'W . (35) 


Br 


ChinaX ive 


222 天 文学 进展 39 3$ 


RRRA MERD) FBT AY Er DLI 7 FR BE SE EIE: 
Sorax = obit B+kr?C gait Agbor B oasr Agbar B 


a, = 0.604 821 875 310 38, a2 = 0.395 178 124 689 62 


36 
k = 0.016 442 830 469 10, b, = 0.224 436 774 742 75 d 
b> = 0.697 313 965 629 25, b3 = 0.078 249 259 628 00 
All 
Sorsy = ptiTAgbir B-Fkr?C qas 7 Agbar B past B 
a, = 0.172 707 170 297 61, az = 0.581 906 819 899 21 (37) 


à3 = 0.245 386 009 803 18, bı = 0.437 551 136 178 33 
bə = 0.562 448 863 821 67, k = 0.011 172 965 730 361 


李 荣 5) mee ee BUS. quxgcs. FAAN, Hénon-Heiles 系统 以 及 量子 力 
学 模型 中 对 新 构造 的 三 阶 最 优化 力 梯度 辛 算法 与 其 他 辛 算法 进行 数值 比较 。 在 经 典 力学 模 
型 中 ， 该 算法 的 能 量 误 差 计 算 占有 绝对 的 精度 优势 ， 并 且 可 以 准确 地 识别 Hénon-Heiles A 
统 的 混沌 轨道 ， 因 此 可 以 有 效 避 免 煞 值 误差 引 起 的 虚假 混沌 现象 。 在 定 态 薛 定 协 方程 的 能 量 
本 征 值 问题 中 ， 该 算法 同样 有 明显 的 精度 优势 ， 甚 至 与 四 阶 标准 辛 算法 相 比 具有 更 高 精度 。 
(2) 四 阶 最 优化 力 梯度 六 算法 
Omelyan 等 人 ™™” 构 造 的 对 称 的 四 阶 最 优化 力 梯度 辛 算法 如 下 : 


OFA: e87B— 0007 Cod rA S 37 B+ allg CeatAgaTB -oT Igi = eW . (38) 
ps ag 7 将 得 到 的 三 阶 最 优化 力 梯度 辛 算法 以 对 称 组 合 的 方式 得 到 两 个 四 阶 辛 算法 ， 


bo a2 a1 
OFAX : Qmd pH oA Tete C gar Ag 7 B4 kr? CsA. TB GTA — erW : (39) 


和 


ag a2 
OFAY : ed rBo BtAg FTE rA bs BA Ir? Ce? TAG BIB, PTA HTB 一 e'W . (40) 


在 开 普 勒 二 体 问 题 和 摄 动 二 体 问 题 等 经 典 力学 模型 ， 以 及 量子 力学 模型 中 的 Morse 势 模型 
中 进行 数值 分 析 比 较 ， 新 的 四 阶 力 梯度 算法 表现 出 了 比较 好 的 数值 精度 ， 甚 至 还 要 明显 优越 
于 已 有 的 四 阶 最 优化 力 梯 度 辛 算法 OF4” 。 

旋转 坐标 系 中 的 平面 圆 型 限制 性 三 体 问 题 的 动能 不 是 动量 的 严格 二 次 型 ， 而 是 受到 旋 
转 坐 标 系 影响 存在 动量 与 坐标 的 交叉 项 ， 这 给 力 梯度 算 子 的 加 入 带 来 了 困难 。 陈 云龙 和 伍 
坎 一 严格 论证 了 力 梯 度 算 子 在 旋转 坐标 系 中 的 形式 仍 与 质心 坐标 中 相同 ， 均 是 引力 的 梯度 
而 不 是 引力 与 非 惯性 力 所 得 合力 的 梯度 ， 只 是 将 A 算 子 变 成 了 如 下 算 子 ， 式 (26) 中 的 C 算 
子 仍然 适用 ， 


4 


-— ð ð ð a a ð 

A p? (n. -T.a.) = (ps t y)g- + (py m + Pug, E - (41) 
应 用 四 阶 力 梯 度 辛 算法 、 最 优化 四 阶 力 梯 度 辛 算法 和 Forest-Ruth 辛 算法 分 别 求解 该 问题 ， 
数值 结果 显示 最 优化 四 阶 力 梯度 辛 算法 能 够 取得 最 好 精度 。 


(ERAT 
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5 隐 式 辛 算法 


5.1 ” 隐 式 辛 算法 的 构造 
当 系 统 的 哈密 顿 变量 不 可 分 离 时 ， 显 式 算 法 一 般 无 法 直接 应 用 。 考 虑 不 可 分 离 哈 密 顿 
函数 : 


H(p, q) = Ho(p, 9) + Hi(p,q) , (42) 
与 3.2 节 中 的 约定 类 似 ， 约 定 Lie SEE A = {, Ho} RH B = {, Hy}. 


正则 方程 表示 为 : pon 
| ge aA = f (p,a) 


p= doo) = g(p,q) = 
ôq 
以 时 间 7 为 步 长 ， 对 整个 哈密 顿 系统 按照 一 阶 Euler BAA Ash, 18580 — BT E o 
辛 算法 M1 和 M*1、 二 阶 隐 式 中 点 法 M2": 
(1) 一 阶 位 置 隐 式 Euler 法 


MI i dm 十 1 = Qn. + TP (Dns qn+1) (44) 
Pn+1 = Pn R5 Tg(Pn; Qn41) 
(2) 一 阶 动量 隐 式 Euler 法 
Mitt. den T da t Tfno Proti) (45) 
Pnt1 = Dn  Tg(dn; Pn41) 
(3) 二 阶 隐 式 中 点 法 
E ; PnTPn+1 Qnin 
Qnt1 =n Tf 2 , 5 
M2: (46) 
p 1=p 十 Tg Pn TPn+1 Int1T In 
ae 2 °” 2 
(4) 二 阶 显 隐 混 合 辛 算法 
A 算 子 是 有 解析 解 的 算 子 ， 而 B 是 隐 式 中 点 法 迭代 求解 的 算 子 ， 两 者 结合 组 成 : 
MSI1 :52 = exp(5A)exp(7B)exp(5B) (47) 
MSI2 :9*2 = exp(7 B) exp(7A) exp(Z B) 


(5) 高 阶 隐 式 辛 算法 
按照 Yoshida" 和 Ruth" 构造 高 阶 显 式 辛 算法 的 思路 对 全 局 哈密 顿 函数 构造 高 阶 隐 式 辛 
算法 : 


| S4(YO) = 52(cr)52(1 — 2cr)S2(cr) (48) 


S4*(YO) = $*2(er)$*2(1 — 2c7)S*2(c7) 


a 
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和 


SA(FR) = A(£r)B(cr)A(45£7)B(1 — 2cr)A(*5£7)B(cr)A($r) 
S4* (FR) = B(§r)A(cr)B(4£r)A(1 — 2cr) B( 
c= 1/2 — 21/8 


4 Fea Nat AE BS Eh RCRA AR, EE Bd AES Ho FAR, Hi 不 可 积 ， 再 根 
HEIR (44) — nere Tower EAA a 算法 的 效率 要 高 于 全 隐 辛 算法 。 

Zhong 等 人 ”将 二 阶 隐 式 中 点 法 及 其 对 称 组 合 运 用 求解 整个 后 牛顿 哈密 顿 系统 ， 数 
值 试 验 表明 ， 显 隐 混 合 中 点 嵌入 法 在 能 量 精度 方面 始终 优 于 全 隐 中 点 法 ， 且 前 者 的 计算 效率 
更 好 。Zhong 和 Wu” 另外 对 中 点 嵌入 法 进行 优化 ， 进 一 步 提 高 了 计算 精度 。 对 于 自 旋 后 午 
顿 哈 密 顿 系统 : 


Hr, p, $1, S2) = Hkep + Hpn . (50) 
开 普 勒 部 分 和 后 牛顿 部 分 分 别 表 示 为 : 
2 
p 1 
Kac e. 
cdd 2 yp (51) 
Hpy = Hipn + Hopn + Hso + Hss 


对 旋转 标量 引入 Wu 和 Xie 的 正则 旋转 坐标 ， 得 到 新 的 正则 哈密 顿 函数 D(r,0,p,£), MEER 
式 中 点 法 直接 应 用 于 求解 变量 D(r0,p,£), WE S2A。 另 外 ，Zhong 和 Wu 分 别 对 开 普 勒 
部 分 和 后 牛顿 部 分 求解 析 解 和 数值 解 ， 数 值 解 部 分 应 用 隐 式 中 点 法 ， 按 照 Yoshida 的 思路 构 
建 四 阶 正则 显 隐 式 混合 辛 算法 : 


S2B = oct odd. 
SAB d. o oif. o odi ue odit 2c)r ZI eie c )7/2 5 off. o og 
A)T 入 Ty 
S4B =p d on ? ma (dk. ep o dL? Q s cdd ? (52) 


(1-23)7/2 。 
ét, o D ES TR © bicep ? 


c =1/2 — 21/53. £ —0.172 086 559 029 5143 , 
x = — 0.091 562 030 755 156 78, A — —0.161 621 762 210 7222 . 


ei 等 人 对 四 阶 Forest-Ruth 型 显 隐 混合 辛 算法 和 四 阶 Yoshida 显 隐 混 合 辛 算法 进 
AU 4 有 二 阶 精度 ， 后 者 可 以 达到 四 阶 精度 。 在 
一 维 不 可 分 系统 也 5(p° +q’) 十 cospsing 和 致密 星 后 牛顿 哈密 顿 系统 中 得 到 了 相同 的 
结论 。 他 们 指出 ， 在 实际 计算 中 ， 当 哈密 顿 的 部 分 变量 可 积 且 解析 解 容易 获得 时 ， 应 采用 
S4(FR) 算法 或 将 解析 解 与 隐 式 数值 解 结合 的 Yoshida 算法 S4(YO)， 以 获得 较 高 的 计算 精度 
和 效率 。 当 可 分 哈密 顿 的 变量 部 分 可 积 ， 但 其 解析 解难 以 求解 ， 或 者 哈密 顿 变 量 部 分 不 可 
积 ， 则 采用 将 解析 解 与 隐 式 数值 解 结合 的 Yoshida 算法 S4(YO). 
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Lubich 等 人 ”提出 了 四 阶 非 正 则 的 显 隐 式 混合 辛 算法 ， 主 要 在 自 旋 致密 双星 后 牛顿 
哈密 顿 系统 中 应 用 ， 将 哈密 顿 分 成 了 三 大 部 分 ， 分 别 是 轨道 Hoa. Ái- Aso FIA 
ie- EE Hss: 


H = Hom + Hso + Hss . (53) 
按照 Suzuki™ 用 二 阶 积分 器 的 五 重 积 构建 四 阶 算法 的 思路 构造 如 下 算法 : 
PR EPH = OES odii 0 GRO” o ORIS O Ori . (54) 


再 分 别 将 这 三 部 分 分 解 为 2, 3, 4 部 分 运用 显 隐 式 混 合 辛 算法 求解 。 将 式 (54) 发 展 到 四 阶 : 
| en = On 9 Onn 9 91, 9 on 9 doa , 
yo = 1/(4— 41⁄8) , y = —41/3/4 — 41⁄3 . 
5.2 ” 几 种 隐 式 辛 算法 的 数值 稳定 性 比较 

当 哈 密 顿 两 部 分 均 可 积 时 ， 刘 福 窑 等 人 ”分 别 对 上 述 一 阶 隐 式 辛 算法 M*1、 隐 式 中 点 
法 M2 以 及 一 阶 显 式 算法 ， 二 阶 显 式 算法 四 种 算法 的 数值 稳定 性 进行 分 析 比 较 ， 将 这 几 种 算 
法 作用 于 线性 哈密 顿 系统 : 

H= yo 上 Do) 1 ETT. Fg) ; U^ z’). (56) 
5 T+V 分 解 相 比 ，Ho + eH 分 解 中 各 算法 的 稳定 区 域 扩大 了 。 2009 年 ， 刘 福 窑 和 钱 晓 
8j"" 比较 了 Liao” 提出 的 两 种 显 隐 混 合 辛 算法 的 数值 稳定 性 ， 即 将 一 阶 隐 式 辛 算法 M*1 和 
隐 式 中 点 法 M2 分 别 嵌 入 到 一 阶 和 二 阶 显 辛 算法 中 得 到 显 隐 混 合 辛 算法 ， 前 者 的 数值 稳定 性 
要 优 于 后 者 。 

钟 双 英 ”分 别 以 一 维 看 合 振子 、 经 典 圆 型 限制 性 三 体 问 题 和 后 牛顿 近似 的 致密 双星 系 
统 为 数值 研究 对 象 ， 对 显 隐 混 合 辛 算 法 MSI 和 MSI2 的 数值 优 劣 性 能 进行 比较 。 一 维 耦 合 
振子 模型 和 圆 型 限制 性 三 体 问题 中 MSI2 算法 的 稳定 性 均 要 明显 优 于 MSI 算法 。 两 种 算法 
在 不 同 的 哈密 顿 分 解 形式 中 表现 有 很 大 不 同 ， 在 荆 十 V 分 解 中 ， 两 种 嵌入 法 都 能 保持 数值 
稳定 ; 在 Ho + Hi DIEP, MSI 算法 在 有 序 轨 道 和 混沌 轨道 中 均 不 能 保持 数值 稳定 。 自 
旋 致 密 双星 后 牛顿 哈密 顿 系统 中 ，MSI1 算法 的 计算 效率 比 MSI2 算法 稍 高 ， 但 稳定 性 不 如 
MSI2 算法 。 综 合 诸 因素 可 知 ， 中 点 舱 入 法 更 适合 于 求解 各 种 相对 论 后 牛顿 动力 学 问题 。 

隐 辛 算法 适用 于 任何 哈密 顿 系统 ， 尤 其 在 使 用 经 典 算法 表现 出 较 差 的 稳定 性 时 ， 隐 辛 算 
法 是 理想 的 选择 ， 但 是 求解 过 程 中 由 于 和 迭代 使 得 计算 效率 大 大 降低 ， 因 此 应 尽量 不 要 作为 首 
选 数值 方法 。 


(55) 


6 扩展 相 空间 显 辛 算法 


Pihajoki^ 提出 一 种 相 空 间 扩充 方法 。 在 不 可 分 离 的 哈密 顿 系 统 中 分 别 增加 一 组 与 原 空 
间 相 同 的 动量 和 坐标 (5, 9)， 得 到 一 个 新 的 可 分 离 的 哈密 顿 量 : 


T (p, p, q, 9) = Hı (,q) + H»(p, q) ý (57) 
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上 式 右边 的 两 部 分 相互 独立 且 可 积 ， 可 以 分 别 进行 解析 求解 ， 由 正则 方程 得 到 : 


~ + OM, " = 1o 

Q=q= Op = {7, Hı}, P=p= ôq = {gq, Hı} , (58) 
oH m aH 

Q=į= ee iH), Pepe-—— (oH). (59) 


显 式 辛 算法 可 以 应 用 ， 初 始 时 刻 (pq). (p.q) 相等 ， 积 分 过 程 中 随 着 时 间 持 续 两 者 变 得 不 
等 ，Pihajoki 将 二 阶 显 辛 算法 即 式 (8)， 改 造 为 : 


Gy (用 ) = Me” /2ehH2/2 M gh Ha /29hHi /2 一 


M2Q(h/2)P(h/2)Q(h/2)P(h/2) Mi P(h/2)Q(h/2)P(h/2)Q(h/2) 


(60) 


Sž (h) a one his /2 yf eh Bi/26 hH2/2 _ 
2Q(h/2)P(h/2)Q(h/2)P(h/2) Mi P(h/2)Q(h/2)P(h/2)Q(h/2) . 
RKF, M, M 表示 两 组 变量 之 间 的 置换 映射 ， 即 : 


(61) 


OM, m 0 0 

au, QM 0 0 
0 0 fw, B M; 
0 0 Bu, fu, 


此 外 ， 需 借助 投影 算 符 W 将 扩展 相 空 间 的 解 返回 到 原始 相 空 间 ， 


w=( aw dw 0 0 iF (63) 
0 0 Bw Bw 


其 中 au, = am tn, au, = õm, In, Bm, = Bur Tn, Py = Bu, Tn, y = awla, dw = dwn, 
© Bw = Bwln, Bw = BwIn J& 8A nx n 的 对 角 和 矩阵 。 原 本 是 辛 对 称 的 ， 但 是 按照 上 述 方法 投 
= 影 到 原 相 空 间 后 就 不 再 保持 辛 结构 。 
, 6.1 ”高 阶 连续 坐标 动量 置换 相 空 间 扩充 显 式 类 辛 算法 

Liu 和 Wu“ 用 偶数 重 积 构建 了 不 同 于 Yoshida 构建 高 阶 辛 算法 思路 的 显 式 高 阶 类 辛 算 
法 ， 用 六 个 偶数 低 阶 算法 构造 高 阶 算法 S4A: 


Sa = S2(Arh)S2(A2h)S2(Agh) Mi x Sx sh) SS Sh) So h) M. 
Sg = Sa(Arh)S4(Azh)S4(Agh) My x S,(h)S4(5h)S4( 5) M: (64) 
Soj42 = = S5j (Arh )S2; (Ash) S5; (Ash) Mi x Sa; (Ash ) Sz; (A2 h) 52; T 


ài = àz = À = 1/2(2 —2V ED), Ag = 1/2 — 2A 


这 种 构建 模式 由 四 部 分 组 成 : 置换 坐标 ， 没 有 置换 的 三 个 偶数 低 阶 蛙 跳 格式 的 积分 ， 置 换 动 
量 ， 没 有 置换 的 三 个 偶数 低 阶 蛙 跳 格 式 的 积分 ” 。 
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在 平面 加 型 限制 性 三 体 问题 、Chin™ 考虑 的 简单 模型 和 无 旋转 致密 双星 的 后 牛顿 哈密 顿 
系统 中 ， 高 阶 连 续 坐 标 动量 置换 相 空间 扩充 显 式 类 辛 算法 在 数值 精度 上 优 于 Yoshida 的 四 阶 
显 式 辛 算法 、Chin” 的 显 式 辛 算法 以 及 显 隐 式 混合 辛 算法 ， 计 算 效 率 也 远 远 高 于 显 隐 式 混合 
辛 算法 。 这 类 显 式 类 辛 算法 也 适合 求解 后 牛顿 哈密 顿 问题 和 旋转 哈密 顿 问 题 等 坐标 和 动量 
不 可 分 离 的 哈密 顿 系 统 。 

6.2 ”中 点 置换 相 空间 扩充 显 式 类 辛 算法 

Liu 和 Wal 提出 的 高 阶 类 辛 算法 需要 进行 两 次 三 重 积 才能 达到 高 精度 效果 ， 且 在 旋转 
致密 双星 后 牛顿 哈密 顿 系统 的 某 些 轨道 积分 中 会 失败 ”。 骆 俊杰 等 人 ”将 连续 坐标 动量 
置换 改 为 一 次 中 点 置换 ( 原 变量 与 它 对 应 扩展 变量 之 间 的 中 点 )， 每 一 步 积 分 都 调整 原 变量 
及 其 对 应 扩展 变量 的 值 为 他 们 的 中 点 值 ， 即 : 


=q. (65) 


再 结合 Yoshida 的 三 重 积 构造 高 阶 类 辛 算法 SAB: 


Soj42 = Mi 2 四 Sa; (Ash) S3; (Agh) $23 (Arh) . (66) 


Sa = Mia & S5 (Ash) SS (Ag h) S4 (A, h) 
Ar = Ag = À = 1/(2 — AVE) | Ms = 1/2 — 2A 


以 不 考虑 引力 耗 散 的 二 阶 后 牛顿 近似 的 自 旋 致密 双星 为 模型 ， 在 周期 轨道 和 混沌 轨道 

中 ， 新 构造 的 中 点 置换 相 空 间 扩充 法 SAB 的 数值 精度 与 隐 式 中 点 法 IM4、 连 续 坐 标 动量 置 
换 相 空间 扩充 法 S4A 相 比 最 高 。 在 混沌 轨道 中 ， 连 续 坐 标 动 量 置换 相 空间 扩充 法 SAA 在 积 
分 一 半 时 间 时 产生 了 剧烈 的 误差 偏 移 。 由 于 Ha, Ha 的 不 同 ， 在 混沌 轨道 中 ， 使 得 S4A 在 置 
换 过 程 中 两 部 分 哈密 顿 函 数 的 误差 不 停 积 累 导 致 失效 。 
中 点 置换 法 可 以 推广 到 非 保 守 的 引力 耗 散 哈密 顿 系 统 中 ， 分 别 以 引力 耗 散 的 2.5PN 后 
牛顿 近似 的 自 旋 致 密 双 星 、 阻 尼 谐 振子 和 受到 Poynying-Robertson 阻力 的 尘埃 粒子 为 模型 ， 
中 点 置换 法 仍然 保持 了 最 佳 的 优越 性 ， 该 方法 的 计算 精度 、 计 算 效 率 以 及 普 适 性 都 优 于 其 
算法 ， 值 得 被 推广 使 用 。 
6.3 ”扩大 相 空 间 FR 最 优化 算法 

吴 亚 林 等 人 结合 优化 思想 和 扩充 相 空间 的 思想 ， 构 造 了 几 种 优化 的 相 空间 扩充 类 辛 算 
法 ， 并 分 别 在 一 维 不 可 分 可 积 系统 、 平 面 圆 型 限制 性 三 体 问题 、 三 阶 后 牛顿 致密 双星 系统 中 
Sees Sut T E T. 

6.3.1 几 种 扩充 相 空间 类 辛 算法 
(1) 连续 坐标 动量 置换 扩充 相 空间 优化 算法 


tS 
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将 FR 算法 和 其 优化 算法 OFR 中 间 算 子 按照 Liu 和 Wa" 构造 的 高 阶 类 辛 算法 进行 拆 
分 ， 分 别 将 (10) RM (16) 式 改造 为 : 


EFR Mae /2eBcheA(l— oh/26B(1-20)b/2 xe eB(1—20)h/26A( -Oh/26BcheAch/2 . (67) 


EOFR —MsePSheAQ-23)h/25Bxhe AA B( -20--6)h/2 Mf, @BU-20c+8))h/2 g^h 


gBXhgAQ-22)h/2,B&h — (68) 


(2) 中 点 置换 扩充 相 空 间 优化 算法 


按照 骆 俊杰 ”提出 的 中 点 置换 法 将 式 (10), (11), (16), (17) 分 别 改造 为 : 
MFR=M®FR. - 
MYO=M@YO. ie 
M cae (71) 
MSU-MGSU. p 
(3) Tao 保 辛 显 式 算法 


Tao 改造 了 Pihajoki 的 方法 ， 将 不 可 分 哈密 顿 系统 改造 为 : 
T (p, Ď, q, q) = Hı (Ď, q) + H»(p, q) + wu Hs(p. p, q, q) . (73) 


ER wH 部 分 是 Tao 人 为 增加 的 约束 哈密 顿 函 数 ，w 是 控制 新 旧 相 空间 变量 之 间 间 隔 的 参 
Jo Hi 和 A 有 解析 解 ， 上 文 第 二 节 中 给 出 的 算 子 A, BAPTA, Ha 部 分 可 积 。 这 
里 再 令 算 子 B 为 : 


e^ = e3Bgh DS $ B 3 (74) 
其 中 DD 算 子 是 Ha 部 分 的 算 子 。 得 到 新 的 二 阶 显 式 辛 算 法 : 
S,(h) = e34e^Be 2A . (75) 


亚 林 等 人 对 FR 算法 以 及 Suziki 算法 用 Tao 的 方法 进行 改造 得 : 
TYO 一 S2(ah)So(1 = ah) S> (ah) ; 
TSU — S5(a,h)S; (a3 h)S3((1 = 2a1 = 2a2)h) S52(a2h) S52(a1h) 


于 算法 中 没有 加 入 置换 因子 ， 算 法 的 辛 结构 将 不 会 被 破坏 ，wH3s 起 到 约束 和 调节 新 旧 空 
间 变 量 的 作用 ，w 的 选取 对 精度 将 产生 很 大 的 影响 。 
6.3.2 算法 的 数值 检验 
分 别 在 一 维 不 可 分 哈密 顿 系 统 、 平 面 圆 型 限制 性 三 体 问题 和 三 阶 后 牛顿 非 自 旋 致 密 双 星 
系统 中 对 以 上 算法 数值 比较 ， 精 度 由 高 到 低 是 MOF R&S MSU >MFRxX MYO > IM4 > 
EOFR ~ TSU > EFR 守 TYO。 总 的 来 说 ， 优 化 后 的 算法 比 未 优化 的 精度 高 ， 各 算法 与 中 
点 置换 结合 比 与 Tao 法 结合 精度 高 ， 能 达到 优化 四 阶 通常 显 辛 算法 的 精度 ， 且 四 阶 隐 式 辛 
算法 的 计算 效率 最 低 。TSU 算法 与 系数 w AK, w= 100 时 精度 达到 四 阶 精度 。 


Sun 
[一 


(76) 


N 


H 


Ti 


C— 
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6.4 “扩大 相 空 间 的 对 数 哈密 顿 显 式 对 称 法 

对 数 哈 密 顿 辛 算 法 在 解决 高 偏心 率 轨道 问题 时 可 以 通过 调节 步 长 避免 数值 失真 ， 构 造 
过 程 较 一 般 的 时 间 变 换 辛 算法 更 简单 。 扩 大 相 空 间 使 得 显 式 辛 算法 可 以 应 用 于 不 可 分 哈密 
顿 系统 ， 且 精度 理想 。Li 和 Wu， "结合 这 两 种 思想 ， 在 Mikkola 等 人 "提出 的 对 数 
哈密 顿 方法 的 基础 上 加 上 一 个 常数 或 者 函数 ， 再 结合 Pihajoki 的 相 空间 扩充 思想 ， 构 造 了 
适用 范围 更 广 的 扩充 相 空 间 的 对 数 哈密 顿 类 辛 算法 ， 它 的 基本 的 哈密 顿 形式 为 ; 


I'(p. p. q, d, po. qo) = In(T | po | C) In(U | C) , (77) 


KH, pp =—-H, p=t, C NETES BAC>0 04 C —0 BI, Æ Mikkola 以 及 苏 湘 
37" 所 构建 和 使 用 的 对 数 哈密 顿 方法 。 时 间 变 换 函 数 g = T+ po+C Ailg=U+C. MAR bi 
的 动能 和 势能 部 分 可 分 时 ， 一 般 的 显 式 辛 算法 可 以 使 用 ;两 者 不 可 分 时 ， 扩 充 相 空间 的 思想 
可 以 被 考虑 。 
6.41 扩大 相 空 间 的 对 数 哈 密 顿 辛 算法 

根据 系统 显 含 时 间 t 和 不 显 含 时 间 t 可 将 哈密 顿 函数 分 别 扩充 为 : 
(1) 不 显 含 时 间 


V(po. p, D. 49,4) = A (po, p, d) + As(po,P.q) . (78) 
(2) 显 含 时 间 


y (po, P, P, Bo, q, q, qo, qo) = Ai (po, qo, P, q) + Na (fo; qo; D; q) j (79) 


Do, do 为 扩充 相 空间 的 物理 量 ， 与 po, qo 对 应 。 
6.4.2 ”数值 检验 

将 四 阶 显 式 辛 算法 FR、 四 阶 隐 式 中 点 法 IS4 (与 6.3.2 节 中 的 IM4 相同 ) 及 三 种 扩充 相 
空间 算法 S4A, S4B, SAC 分 别 应 用 求解 开 普 勒 摄 动 二 体 问 题 、 圆 型 限制 性 三 体 问题 和 无 自 旋 
三 阶 后 牛顿 致密 双星 系统 。 在 三 种 模型 中 ，S4C 算法 的 表现 始终 是 最 好 的 ， 隐 式 中 点 法 IS4 
算法 虽然 也 表现 出 了 很 好 的 计算 精度 ， 但 是 在 计算 成 本 上 最 高 ，S4B 算法 与 IS4 算法 精度 同 
量 级 且 计 算 成 本 最 低 。 扩 大 相 空 间 的 对 数 哈密 顿 显 式 类 辛 算法 在 数值 稳定 性 、 数 值 精度 以 及 
计算 效率 上 均 占 有 很 大 的 优势 ， 且 适用 于 高 偏心 率 轨 道 。 本 文 主要 在 圆 型 限制 性 三 体 问题 
中 对 以 上 算法 进行 数值 比较 ， 这 几 种 算法 分 别 是 四 阶 显 式 辛 算法 FR、 四 阶 隐 式 中 点 法 ISA 
连续 坐标 动量 置换 相 空间 扩充 显 式 算法 EFR 及 优化 算法 BOFR、 中 点 置换 相 空 间 扩 充 显 式 
算法 MFR 及 优化 算法 MOFR。 如 图 四 所 示 ， 对 于 低 偏心 率 轨道 ， 四 阶 显 式 辛 算法 FR 在 长 
时 间 积分 过 程 中 能 保持 能 量 误差 稳定 且 精 度 较 高 ， 但 在 高 偏心 率 轨道 中 ，FR 算法 在 长 时 间 
积分 后 能 量 误差 开始 增长 。 c=0 时 ， 连 续 坐 标 动量 置换 EFR 四 阶 算法 与 其 优化 算法 EOFR 
均 不 适用 于 高 偏心 率 轨道 ， 它 们 在 较 短 时 间 误 差 就 开始 增长 。c=1.7 时 ， 各 相 空 间 扩 大 法 的 
精度 均 有 提高 ， 且 EFR 算法 和 EOFR ,算法 在 高 偏心 率 轨道 中 使 得 能 量 误差 不 再 线性 增长 ， 
说 明 系 数 c 的 选取 很 重要 。 另 外 ， 中 点 置换 相 空间 扩大 法 MFR 及 其 优化 算法 MOFR 的 计 
算 效率 最 高 ， 四 阶 隐 式 算法 IS4 的 计算 效率 最 低 。 
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ik: 图 中 黑 线 表示 FR. “ZARA EFR, HARA MFR, MARA EOFR. KARA IS4， 黄 线 表示 


图 3 圆 型 限制 性 三 体 问题 中 的 两 条 轨道 的 能 量 误差 


7 总 结 与 展望 


辛 算法 自从 被 提出 以 来 ， 特 别 是 可 分 哈密 顿 系统 的 显 式 辛 算法 ， 在 动力 学 天 文中 被 广 
泛 应 用 ， 已 充分 显示 出 传统 算法 不 可 比拟 的 优点 。 在 定量 天 文学 计算 问题 中 ， 辛 算法 可 


以 有 效 控制 轨道 治 迹 误 差 的 持续 快速 增长 ， 但 是 对 于 同 阶 算法 ， 辛 算法 的 截断 误差 要 比 
Runge-Kutta 法 等 传统 方法 大 ， 在 较 短 时 间 内 优点 无 法 突出 。 辛 校正 和 类 高 阶 辛 算法 则 通过 


误差 校正 ， 在 不 缩小 积分 步 长 的 前 提 下 使 得 算法 精度 提高 ， 时 间 变 化 辛 算法 使 用 变 步 长 ， 解 


决 了 天 体 紧 密 交 会 和 大 偏心 率 轨道 中 的 数值 解 失真 问题 。 
它 也 能 保证 系统 的 辛 结构 和 能 量 误差 稳定 ， 但 积分 过 


算法 非常 适合 应 用 ， 与 传统 算法 相 比 ， 


对 于 不 可 分 哈密 顿 系统 ， 隐 式 辛 


程 由 于 多 次 迭代 使 得 计算 效率 显著 降低 。 扩 大 相 空 间 的 类 辛 算法 近年 来 得 到 了 快速 的 发 展 ， 


其 中 Luo 5 A 7 构造 的 中 点 置换 扩充 相 空间 方法 是 目前 最 理想 的 ， 它 具 


高 精度 、 高 效 
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用 范围 ， 之 后 提出 的 扩大 相 空 间 FR 最 优化 算法 ”以 及 扩大 相 空 间 对 数 哈密 顿 
四周 都 对 此 进行 了 证 明 。 力 梯度 辛 算法 作为 一 种 独特 的 辛 算法 ， 在 时 间 可 首 


的 哈密 顿 系 统 求解 中 具有 很 大 的 优势 ， 但 在 哈密 顿 变量 不 可 分 时 应 用 有 


ES XE. 


参考 文献 : 


o o) Noan AUNDH 


e U Q0 00 Q2 Q2 Q2 Q0 C2 CO Q2 L2 LO Lb PR P P2 P2 M bl 2 Be bhe he hhe FPF RP RP RP RP h 
O Oo AND Te wWwnr OM WAND TF Wwnr OV WAN DWT ek Wwnr O 


Feng K. difference schemes and symplectic geometry. Beijing: Science Press. 1984: 42 


Ruth R D. IEEE Transactions on Nuclear Science, 1983, NS-30(4): 2669 
Forest, Ruth R D. PhyD, 1990, 43(1): 105 

Kinoshita H, Yoshida H, Nakai H. CeMDA, 1991, 50: 59 
Yoshida H. PhLA, 1990, 150: 262 

Wisdom J, Holman M. AJ, 1991, 102: 1528 

Chamber J E, Murison M A. AJ, 2000, 119(1): 425 

Wisdom J, Holman M, Touma J. Fields Institute Communication, 1996, 10: 217 
Chin S A. PhLA, 1997, 226(6): 344 

Chin S A, Chin C R. Celest. Mech Dyn. Astron, 2005, 91: 301 
Mclachlan R I, Atela P. Nonlinarity, 1992, 5: 541 

ER. 硕士 论文 . 南昌 : 南昌 大 学 , 2011 

ÆR, (EK. 物理 学 报 , 2010, 59(10): 7135 

Li R, Wu X. SCPMA, 2010, 53(9): 1600 

陈云 龙 , thik. 物理 学 报 , 2013, 62(14): 41 

Liao X H. CeMDA, 1997, 66: 243 

Suzuki M. PhLA, 1990, 146: 319 

Wu X, Xie Yi. PhRD, 2010, 81: 084 

Zhong S Y, Wu X. PhRD, 2010, 81: 1 

Zhong S Y, Wu X, Liu S Q, et al. PhRD, 2010, 82: 124 

Wu X, Zhong S Y. GeReGr, 2011, 43: 2185 

Mei L J, Ju M, Wu X, et al. MNRAS, 2013, 435(3): 2246 

Mei L J, Mei L J, Liu F Y. EPJC, 2013, 73(5): 2413 
Pihajoki. CeMDA, 2015, 121 (3): 211 

Tao M. PhRE, 2016, 94: 043 

Liu L, Wu X. MNRAS, 2016, 459(2): 1968 

Liu L, Wu X. Gen Relativ gravit, 2017, 27(4): 397 

Luo J J, Wu X, Huang G Q, et al. ApJ, 2017, 834: 64 

骆 俊杰 . 硕士 论文 . 南昌 : 南昌 大 学 , 2017 

Wu Y L, Wu X. IJMPC, 2018, 29(1): 6 

Li D, Wu X. MNRAS, 2017, 469: 30 

Mikkola S, Palmer P, et al. CeMDA, 1999, 74: 289 

Preto M, Tremaine S. AJ, 1999, 118: 2532 

Omelyan I P, Mryglod I M, Folk R. CoPhC. 2002, 146(2): 188 
Wisdom J, Holman M. AJ, 1992, 104: 2022 


Varadarajan V S, Lie group. Lie algebras and their representation, Englewood Cliffs: Perntiee Hall, 1974 


Blanes S, Casas F, Murua A. MNA, 2008, arXiv: 0812 
Deng C, Wu X, Liang E W. MNRAS, 2020, 496: 2946 
Mikkola S, Palmer P. CeMDA, 2000, 77: 305 

Duncan M J, Levison H F, Lee M H. AJ, 1998, 116: 2067 


刊 


232 天 文学 进展 39 3$ 


41] Levision H F, Duncan M J, Zahnle K, et al. AJ, 2000, 120: 2117 
42] Malhotra R. Icarus, 1991, 94: 399 

43] Wu X, Huang T Y, Hong Z, et al. Ap&SS, 2002, 283: 53 

44| Wu X, Huang T Y, Wan X S. A&A, 2003, 27: 114 

45] Quinlan G D, Tremaine S. AJ, 1990, 100: 1964 

46] 黄 天 衣 , EBH, 赵 长 印 . 天 文学 报 , 1997(03): 278 

47| Laskar J, Robutel P. CeMDA, 2001, 39: 52 

48] Chin S A, Chin S R. Chem Phys, 2001, 114(17): 7338 

49] Chin S A, Anisimov E. Chem Phys, 2006, 124(5): 054 

50] Chin S A. PhRE. 2007, 75: 036 

51] Sun W, Wu X, Huang G Q. RAA, 2011, 11(3): 353 

52] Xu J, Wu X. A&A, 2009, 10(2): 173 

53] 徐 佳 . 硕士 论文 . 南昌 : 南昌 大 学 , 2010 

54] Omelyan I P, Mryglod I M, Folk R. PhRE, 2002, 66(2): 026 
55] Omelyan I P, Mryglod I M, Folk R. CPhC, 2003, 151: 272 
56] Li R, Wu X. EPJP, 2011, 126: 73 

57] FR. 安阳 师范 学 院 报 , 2015, 5: 6 

58] Lubich C, Walther B, Brigmann B. PhRD, 2010, 81(10): 1292 
59] Xia, fl A, 陆 本 魁 . 天 文学 报 , 2006, 4: 418 

60] 刘 福 窗 ,， 钱 晓 明 . 南昌 大 学 学 报 (理科 版 ), 2009, 33(02): 123 

61] 钟 双 英 . 硕士 论文 . 南 南昌 大 学 , 2011 

62] RER, 硕士 论文 . 南昌 : 南昌 大 学 , 2018 

63] 李 丹 . 硕士 论文 . 南昌 : 南昌 大 学 , 2018 

64] Mikkola S, Tanikawa K. New Astronomy, 2013, 20: 38 

65] Mikkola S, Tanikawa K. MNRAS, 2013, 430: 2822 

66] HMT. 硕士 论文 . 南昌 : 南昌 大 学 , 2016 


Iu Inu 


Classification and Development of Symplectic Algorithm 


SUN Lang!?, LIU Fu-yao!, WANG Ying!?, SUN Weil? 


(1. Shanghai University of Engineering Sciences, School of Mathematics, Physics and Statistics, Shang- 
hai 201620, China; 2. Shanghai University Of Engineering Sciences, Center of Application and Research 
of Computational Physics, Shanghai 201620, China; 3. Anqing Normal University, School of Resources 
and Environment, Anqing 246133, China) 


Abstract: Symplectic algorithms are the best integration tool to study the long-term evolu- 
tion of Hamiltonian systems. Because of this good property, they have been popular. One can 
construct high-precision symplectic algorithms by analyzing the truncation error of Hamilto- 
nian function. Symplectic algorithms consist of explicit symplectic algorithms and implicit 
symplectic algorithms. When a Hamiltonian system is integrable or separable, explicit sym- 
plectic algorithms can be constructed. Explicit algorithms includes nonforce gradient explicit 


symplectic algorithms, force gradient symplectic algorithms, symplectic correction methods 
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and pseudo high order symplectic algorithms. When the variables of a Hamiltonian system 
cannot be separated, implicit symplectic algorithms, or the explicit symmetry methods in 
extended phase space are suitable for solving this problem. We introduce the construction 
forms and applications of these algorithms, analyze the advantages and disadvantages of the 
related symplectic algorithms. It is suggested that symplectic algorithms should be selected 


according to practical problems. 


Key words: Symplectic algorithm; Hamiltonian system; Symplectic correction; force gradi- 
ent symplectic algorithm; explicit symplectic-like algorithm in extended phase 


space 


